A molecular dynamics "Maxwell Demon" experiment for granular mixtures* 
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We report a series of molecular dynamics simulations and investigate the possibility to separate a 
granular mixture of inelastic hard spheres by vigorously shaking it in a box made of two connected 
compartments. As its one-component counterpart, the system exhibits a "left-right" symmetry 
breaking entirely due to the inelasticity of grain-grain collisions, and triggered by increasing the 
number of particles. In the compartment where the density of grains is larger, we observe a partial 
segregation with a predominance of heavy particles. However, this compartment still has a higher 
density of light particles than the other one, which is light-rich. The density, granular temperature 
and anisotropic pressure profiles are monitored. We also discuss how to construct a relevant order 
parameter for this transition and show that the resulting bifurcation diagram is dominated by large 
fluctuations. 



I. INTRODUCTION 



C3 



-a 
c 

o 
o 



> 
o 

(N 

(N 
O 

-i— > 
c3 



i 

G 

O 

o 



X 



Although granular matter may exhibit similarities with molecular fluids (such as pattern formation), it is neverthe- 
less intrinsically out of equilibrium: The inter-particle collisions dissipate kinetic energy, and a steady state may only 
be achieved by a suitable energy supply. As a result, such systems may display many phenomena that are "forbidden" 
by the laws of equilibrium statistical mechanics. In the realm of granular gases (dilute systems of macroscopic grains 



|4], non Gaussian velocity distribu- 
L6ril7], and breakdown of kinetic energy 



in rapid motion and colliding inelastically), tendency to form clusters 
tions §, |, @, |, [ll], H || , long range velocity correlations |J| || _ 

equipartition in a mixture of dissimilar grains [l^, 19, 20 , |TJ have been reported. 

Another interesting feature, at complete variance with equilibrium phenomenology has been obtained with a simple 
experiment |2^, ^4|, : a vibrated system of grains confined in a box with two connected identical compartments 
may exhibit a stationary state with spontaneous symmetry breaking (non-equipartition of grains between the two 
compartments). This clustering phenomenon may be interpreted as a separation in a "hot" and a "cold" region, 
considering that the granular temperature is a direct measure of the mean squared velocity of the particles. In 
the limit where the exchange of particles between the two compartments may be considered as an effusion process, 
Eggers [^6| has put forward an analytical approach to explain this apparent intrusion of a "Maxwell Demon". On 
the other hand, Brey and collaborators 27 reported a hydrodynamic mechanism for the symmetry breaking, that 
becomes operational under some simplifying assumptions in the opposite limit where the size of opening connecting 
the two compartments is larger than the mean-free-path of the gas in its vicinity. 

In this contribution, we revisit numerically the Maxwell Demon experiment in the latter case, and consider the 
specific situation of a binary low density granular mixture, with the aim to investigate whether such a set-up is able to 
achieve an efficient segregation of t he m ixture. The model is defined in section Making use of molecular dynamics 
simulations, we discuss in section III how to construct a relevant order parameter for the transition under study, 
and show that it is dominated by large fluctuations (as also observed recently in a related context 28 1). The two 
components of the mixture are found to behave differently: heavy particles display a stronger (left-right) asymmetry 
than the light ones, leading to a separation between a dense gas rich in heavy particles and a dilute light-rich gas. 
The behaviour of partial densities and granular temperatures are investigated (section IV ) from which we deduce the 
different components of the pressure tensor making use of the general equation of state derived in jl3) ; although the 
hypothesis of an isotropic pressure given by the ideal gas equation of state is clearly not verified, we show that the 
no-convection hydrodynamic condition of a divergence-free pressure tensor (V • P = 0) is obeyed, taking into account 
anisotropies, boundaries and corrections to the ideal gas equation of state. Conclusions are finally drawn in section 
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II. THE MODEL 

The system is made of N inelastic hard disks evolving in a 5 x L two-dimensional box, losing energy at inter- 
particle collisions and gaining energy through collisions with two vibrating walls situated at y = and y = L. The 
particles have diameters er^ and masses rrij, i = 1, 2. A binary collision between grains of species i and j is momentum 
conserving and dissipates kinetic energy: the collision i-j is characterized by the coefficient of normal restitution a^j. 
Accordingly, the pre-collisional velocities (Vi,Vj) are transformed into the post-collisional couple (v^v'j) such that 

v\ = vi -f- — (1 + Oij){a ■ Vij)a (1) 

rtii + rrij 

where v t j — Vi — Vj and a is the center to center unit vector from particle i to j. Note that cty = au to ensure 
the conservation of total linear momentum rriiVi + rrijVj. The total density is denoted p, and the partial densities 
Pi = %%p (the number fractions Xi are such that J2i x i = !)■ The granular temperature of species i is Tj, defined from 
the mean kinetic energy of subpopulation j, by analogy with the usual temperature of elastic gases: Tj = (rriivf)/d, 
where d is the space dimension (here d = 2). In the remaining of the paper this granular temperature will be coined 
"temperature" for simplicity. 

The box is divided into two compartments of width 5/2 by a wall parallel to Oy starting at height y$. The walls 
located at x = and x = 5 are elastic, while those at y — and y = L are vibrating and thus inject energy into the 
system. For simplicity, the two vibrating walls are taken to move in a saw-tooth manner, so that a colliding particle 
at y — (resp. y = L) always finds the wall to move "upwards" (resp. "downwards") with the same velocity v$ (resp 
—vq). In addition, the amplitude of the vibration is considered to vanish (i.e. to be much smaller than the local mean 
free path ^7j), so that the walls are located at the fixed positions y — and y = L: the y-component velocity 
of a particle colliding with the wall at y — (resp. y — L) is therefore changed according to v' y = 2v n — v y (resp. 
v' y = —2vq — Vy). Since we consider vigorous shakings, the gravitational field has not been included in the analysis. 

For simplicity, we have considered equimolar mixtures (N\ = N2) of particles having the same diameter (cti =172) 
but different masses. Various mass ratios mi/m2 S [1 : 5] have been studied |^|, so that the species 1 is always the 
heavier particle. We have run molecular dynamics simulations changing N either at constant packing fraction 
(equal to 7r j O(T 2 /4 in two dimensions) or at constant <Tj, with the same qualitative observations. The numerical results 
we will present correspond to a fixed low mean packing fraction t]q = 0.015 (the inelastic collapse H occurring if 
the mean density exceeds a low threshold), and to equal coefficients of restitution ctij = 0.9, close to experimentally 
relevant values. We have investigated other values of the restitution coefficients between 0.7 and 0.9, and two different 
aspect ratios, L = S and L — 25, with the same qualitative results. 

III. BIFURCATION DIAGRAM AND LARGE FLUCTUATIONS 

For a one component system, it has been shown from a hydrodynamic approach [ p7[ that, as the number of particles 
in the box (N) is increased, a transition occurs at a certain threshold N*: for N < N*, the system is symmetric, i.e 
the mean number of particles in each compartment is AV2 while, for N > N*, one of the compartments becomes more 
populated and colder than the other. This hydrodynamic study relies of the assumption of an isotropic pressure given 



by the ideal gas equation of state 27 . At a given inelasticity (i.e. at given values of the restitution coefficients), the 
control parameter (governing the transition from symmetric to asymmetric situation) is proportional to Na^ 1 j S , 
where a is the particle diameter [ p7[ . At fixed reduced density n — Na d /(LS) the above parameter is proportional 
to N x / d , while at fixed size a it scales like N. 

The "order parameter" of this transition was defined in |2^, |27j as the time average (e) of the asymmetry e: 

N-2N le f t 

e = 3 

27V V ; 

where N le ^ 1 is the number of particles in the left compartment. In addition to the global e, we may introduce two 
relevant asymmetry parameters for each type of particles: 

N- — 2N [eft 

— < = 1 ' 2 - (4) 

For a given simulation time, if one computes (cj) for the binary mixture (or |(ej)| to have a positive quantity), a 
left-right symmetry breakdown is evidenced (see Fig. 0). The asymmetry is more pronounced for heavy particles 
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(|(ei)| > | (e2 ) | ) , and |(ei)| increases with the mass ratio mi/rri2- On the other hand, the light particles asymmetry 
decreases with mi/m 2 . At this point we conclude that the compartment with larger global density is heavy-rich, 
while the lighter particles are more uniformly distributed and therefore the less populated compartment is richer in 
light particles. 

However, for symmetry reasons, one should expect that the mean value (e) (and the (e^)) always vanish for sufficiently 
long simulation times, so that these quantities are arbitrary and do not provide relevant order parameters. Inspection 
of the time behaviour of N^^ 1 confirms this picture [see Fig. ^-a)], which is made more quantitative by computing 
the probability distribution function of e over very long runs and various initial conditions [see Fig. g-b)]. At small 
N, e fluctuates around and its standard deviation increases with N. As N increases, the asymmetric configurations 
become stable but the system continuously jumps from one of the possible asymmetric situations to the other, still 
spending some time in between close to the symmetric state. When N further increases, the residence time spent 
in each of the asymmetric states increases and may eventually overcome the simulation time: for N 3> N*, starting 
from a symmetric situation, the system quickly evolves into an asymmetric configuration, in which one compartment 
is strongly overpopulated, and remains in this situation for all the simulation time. For larger simulation times 
however the symmetry would be restored. The situation is thus analogous to that of a two-states system, in which 
the energy barrier between two symmetric states increases with system size. This behaviour is reminiscent of that 
recently reported in [psj: in this study of a translational symmetry breaking as the aspect ratio of the simulation 
box is changed (without a separating wall), large fluctuations have been shown to occur in a wide region around a 
hydrodynamically predicted threshold value beyond which the homogeneous system becomes unstable. 

In any case, | (e) | vanishes for any N for long enough simulations, and does not provide an acceptable order parameter. 
There are then a priori two possibilities to construct such a quantity: a) by time averaging |e| or b) by extracting the 
most probable (say positive) value e* of e from its probability distribution function (p.d.f, see Fig. ^|), averaged over 
the simulation time and over various initial conditions. Note that this p.d.f., and thus both possible definitions, are 
not sensitive to the length of the simulations (except for very small simulation times) |3lJ]. We compare in Fig. g these 
two definitions with the previous one, |(e)|, computed again for a given (large) simulation time. Since e fluctuates 
around even at small N, the (|ej|) depend rather smoothly on TV, and therefore do not allow a clear definition of a 
critical number of particles. On the other hand, the most probable values e* allow to define a critical region, being 
identically at small N and taking positive values above a certain threshold (see Fig. g). When N is large enough, 
the probability distribution functions of the become sharply peaked around the e* so that both quantities e* and 
become close; moreover, these p.d.f. take extremely small values in the vicinity of ei = 0: this corresponds to 
the fact that the system is stuck for long times in one of its two most probable states, so that the computation of 
|(ej)| coincides with that of (|ej|), or e*. 



IV. DENSITY PROFILES AND PRESSURE TENSOR 

Instantaneous typical configurations are displayed in Fig. ^| for various values of N and mass ratio m\jmi- From 
the coarse grained local packing fractions 7]i(x,y), we define x-averaged quantities in each compartment: 

2 f s / 2 2 f s 

n l i{y) = -5! dxrn(x,y) , rfi(y) = - / dx rn(x,y) ,i=l,2. (5) 

JO ° JS/2 



These quantities are averaged over time for one run between two successive "flips" (see section [IIlD , but not averaged 
over various runs since the asymmetry would then be lost. The corresponding density and temperature profiles are 
shown in Fig. g for N = 1000, well above the bifurcation point. One may observe that in the asymmetric situation, 
the densities are different even for y < yo, i.e. not only where the compartments are physically separated. 

Two-dimensional plots of the coarse grained densities f]i(x, y) are displayed in Fig. |?j for two values of the number 
of particles, well below and well above the bifurcation. Below the transition, translational invariance in x holds in the 
whole box, while above, the densities and temperatures are almost independent of x in each compartment, but are 
discontinuous at x = 5/2 for y > yo because of the separating wall; at y < yo but close to yo a quite sharp change is 
observed in the vicinity of x — S/2. At small y the x-gradients are smaller. 

In the hydrodynamic study of j27j, the ideal gas form for the pressure constitutes a fundamental hypothesis which 
allows for an analytic treatment; moreover, the pressure is assumed to be isotropic. However, anisotropic energy 
injection mechanisms lead to anisotropic pressure tensors, especially near vibrating walls [^, [32|. Knowing the 
density and pressure profiles for our system, one may compute the two components P xx and P yy of the pressure tensor 
from a given equation of state. We consider the generic expression derived in |]l3|| within Enskog-Boltzmann kinetic 
theory. For an homogeneous and isotropic mixture with partial temperatures Ti, number fractions Xi, and without 
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any approximation on the single particle velocity distribution, it was obtained that 

^ d 

t— 1 ^— ' mi + rrii (cr ) 

where d denotes space dimension, cry = (<jj + <Tj)/2, (a d ) = ^Ylii x i. <T i- The ^ij are the a priori unknown pair 
distribution functions at contact; these quantities embody the correction to the ideal gas equation of state, and since 
we are considering dilute system, it is sufficiently accurate to assume Xij — 1 (l° w density limiting value). The values 
of P xx and P yv are finally obtained by substituting T,- in @ respectively by Ti X and Ti V (i.e. the mean square x or y 
components of particles velocities) . 

The results are summarized in Fig. |^, where we plot P xx (x,y) and P yy {x,y) for N — 600 and ./V = 900. Below 
the transition (TV = 600), the picture is similar to the one without a separating wall, and the pressure tensor is x 
independent. One therefore has d x P xx = d x P yy = d y P yy = 0. However, for TV = 900 (above the transition), the yy 
components are no longer equal in the left and right sides, while the xx components are equal only for y < yo where 
the separation begins. One still observes d x P xx = d y P yy = in each compartment, except close to the extremity of 
the separating wall (x — S/2,y — yo), but for y > yo the separating wall allows for different values of the pressure 
components. Moreover d x P yy ^ for y < yo while for y > yo, d x P yy — in each compartment (except close to 
x = S/2), with a discontinuity at x = S/2. It is noteworthy that the denser compartment is also the one where both 
components of the pressure tensor are lower, since it is much "colder" than the dilute compartment. 

This analysis shows that both above and below the symmetry breaking, the (anisotropic) pressure tensor as com- 
puted from Eq. (|^) is divergence free: d x P xx +d y P yy — 0. This "hydrostatic" requirement follows from the condition 
of a vanishing flow field, and in spite of the low mean densities considered here, would not be fulfilled restricting P xx 
and P yy to their ideal parts. 



V. CONCLUSIONS 



For a one component granular gas enclosed in a box made of two connected compartments, a vigorous shaking is 
known to promote a symmetry breakdown and separate the system into a cold and dense region on the one hand, and 
a hot and dilute part on the other hand (so called "Maxwell-Demon" experiment). In addition, in a binary granular 
mixture, heavy and light grains generically have different granular temperatures. In this contribution, we have 
combined both aspects (Maxwell Demon and mixture) to investigate the possibility to separate the two components 
of the mixture. Our molecular dynamics results show a spontaneous symmetry breaking as the number of particles 
is increased, all other parameters being kept constant. The denser compartment then appears to be rich in heavy 
particles, but this partial segregation is such that this compartment is also richer in light particles than the other half 
of the confining box (which is however the light-rich one). It therefore seems that such a set-up cannot achieve an 
efficient segregation (although a possibility would be to isolate the dense compartment and iterate the process with 
this non-equimolar mixture). 

The transition reported here is not stricto sensu a phase transition, since the control parameter is the system size. 
As a consequence, fluctuations can always bring the system from one of the asymmetric states to the other, as e.g. 
for a finite size Ising model below its critical temperature. We have discussed the consequences of this feature on the 
definition of a relevant order parameter to characterize the bifurcation. 

We finally note that the experimental realization of the two dimensional situation investigated here seems feasible, 
for instance by adapting the configuration used in |(| (friction with the walls confining the system in a 2D slab might 
play an role, and has not been considered here) . The experimental signature of the large fluctuations that invalidate 
hydrodynamic approaches seems an interesting issue. 
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FIG. 1: Schematic picture of the set-up 
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FIG. 2: (a): Asymmetry parameters |(e)| and |(e»)| versus number of particles for three different mass ratios and a given 
simulation time, i.e. number of collisions per particle. Here, all the inelasticity parameters are taken equal: an = ai2 = 022 = 
0.9. The opening connecting both compartments is 40% of the total height of the simulation cell (yo = OAL). 
(b): Asymmetry parameters versus mass ratio, at fixed number of particles N = 1000 and a%j = 0.9. As the mass ratio 
increases the asymmetry increases for the heavy particles and decreases for the light ones. As emphasized in the text, these 
figures depend on the simulation time available. 
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FIG. 3: (a) Number of particles of type 2 in the left compartment as a function of time (measured in number of collisions per 
particle), for various values of N, and mi/m.2 = 3. In all cases, the horizontal lines correspond to the symmetric situation 
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(b) Probability distribution function of e for rri2/mi = 3. 
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FIG. 4: (a) Asymmetry parameters (|ei|) and |(ei)| versus number of particles for mi = 5m2. The inelasticity coefficients are 
the same as in Fig. ^ (ptij = 0.9). 

(b) Comparison of the most probable values e* with (|ei|) for the same parameters as in (a). 
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FIG. 5: Typical instantaneous snapshots. Heavy particles (label 1) are denoted by a plus, and light particles (label 2) by a 
circle. Top left: N = 650, mi = 3m,2; Top right: N = 850, mi = 3ni2; Bottom left N = 1000, mi = 2m,2; Bottom right: 
AT = 1000, mi = 5m 2 . 
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FIG. 6: (a) Density and (b) temperature profiles for N = 1000, mi = 2ui2- The left compartment (0 < x < S/2) is denser and 
colder than the right one. In the right compartment, the light particles are denser than the heavy ones. The mean packing 
fraction, averaged over the whole system is t)q = 0.015. The ratio r)(y)/ijo is also the ratio p{y)/p of local density normalized 
by the mean one. The separation between the two compartments is located at x = S/2; 0.4 < y/L < 1. 




FIG. 7: Averaged local density of the particles of type 2 (light component) for N = 400 (left panel) and N = 900 (right panel). 
The separating wall is at x — S/2, y > y = 0AL and mi = 3m,2 in both cases. 
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FIG. 8: Components of the pressure tensor P xx (left) and P yy (right), as given by the equation of state [Eq. (g)], for N = 600 
(top) and N — 900 (bottom). Here, mi/ma = 3 and cuj ~ 0.9. For N = 900, the right compartment (a; > S/2) is more 
populated, colder and at a lower pressure. 



